package modeling;


/**
 * @author Mathieu Rousseau
 * @email mathieu.rousseau3@mail.mcgill.ca
 */


import java.io.BufferedReader;
import java.io.BufferedWriter;
import java.io.File;
import java.io.FileReader;
import java.io.FileWriter;
import java.io.IOException;
import java.text.DateFormat;
import java.text.DecimalFormat;
import java.text.NumberFormat;
import java.text.SimpleDateFormat;
import java.util.Date;
import java.util.Random;

import utilities.PDButility;


/** ------------------------------------------------------------------------------------------------------------------
 * start of MCMC5C class
 * ------------------------------------------------------------------------------------------------------------------  **/
public class MCMC5C {


	/*  global variables  */
	private Date date;
	private DateFormat dateFormat = new SimpleDateFormat("yyyy/MM/dd\tHH:mm:ss");
	private final DateFormat dateFormatOutputFolder = new SimpleDateFormat("yyyy-MM-dd-HH-mm-ss");
	private final String outputFolderDate = dateFormatOutputFolder.format(new Date());
	private final String outputFolder = "MCMC5C/Output/" + outputFolderDate;
	private final String inputPath = "MCMC5C/Input/";
	private final String outputPath = "MCMC5C/Output/" + outputFolderDate + "/";
	private NumberFormat formatter = new DecimalFormat("0");
	private PDButility pdb = new PDButility();
	private String filenameIFs = "no_name.txt";
	private String filenameFragments = "no_name.txt";
	private String initialStructureFilename = "no_name.txt";
	private int initialStructureCase = 1;
	private int numberOfPoints = 0;
	private boolean determinedNumberOfPoints = false;
	private int[] pointGenomicLocations;
	private int numberOfIFValues = 0;
	private int numberOfVarianceValues = 0;
	private double[][] IFvalues;
	private double[][] IFvariances;
	private double averageIF = 0;
	private double[][] structure;
	private double[][] bestStructure;
	private double probabilityOfBestStructure;
	private long indexOfBestStructure;
	private long runLength = (long)1E7;
	private double numberOfStructuresToOutput = 10;
	private double stepSize = 0.05;
	private double exponent = 1;
	private long numberOfChangesAccepted = 0;
	private String[] argsGlobal;
	/*  end of global variables  */
	
	
	/*  start of getNumberOfPoints method  */
	private void getNumberOfPoints(String fileToCountNumberOfPointsFrom){
		String line = null;
		numberOfPoints = 0;
		try{
			BufferedReader br = new BufferedReader(new FileReader(inputPath + fileToCountNumberOfPointsFrom));			
			while((line = br.readLine()) != null){
				numberOfPoints++;
			}
			br.close();
		}catch (IOException e) {
			e.printStackTrace();
			System.out.println("Error: cannot find file with filename \"" + fileToCountNumberOfPointsFrom + "\"\n" +
			"System.exit(0)");
			System.exit(0);
		}
		determinedNumberOfPoints = true;
	}/*  end of getNumberOfPoints method  */

	
	/*  start of getFragmentsFromFile method  */
	private void getFragmentsFromFile(){
		//  make sure that the number of points has already been determined
		if(!determinedNumberOfPoints){
			System.out.println("Error: the number of points must be determined before reading the fragments as input.\nSystem.exit(0)");
			System.exit(0);
		}
		
		pointGenomicLocations = new int[numberOfPoints];

		String line = null;
		int lineCounter = 0;
		try{
			BufferedReader br = new BufferedReader(new FileReader(inputPath + filenameFragments));			
			while((line = br.readLine()) != null){
				String values[] = line.split("\t");
				if(values.length == 1){
					pointGenomicLocations[lineCounter] = Integer.parseInt(values[0]);
				}
				else{
					System.out.println("Error: there is not one entry in each row in the fragments file\n"
							+ "System.exit(0)");
					System.exit(0);
				}
				lineCounter++;
			}
			br.close();
		}catch (IOException e) {
			e.printStackTrace();
			System.out.println("Error: cannot find fragment file with filename \"" + filenameFragments + "\"\n" +
			"System.exit(0)");
			System.exit(0);
		}
		if(lineCounter != numberOfPoints){
			System.out.println("Error: there are not the same number of fragments (" + lineCounter + ") as there are points (" + numberOfPoints + ")\n"
					+ "System.exit(0)");
			System.exit(0);
		}
		System.out.println("Opened and read file: " + filenameFragments);
	}/*  end of getFragmentsFromFile method  */


	/*  start of getIFsFromFile method  */
	private void getIFsFromFile(){
		//  make sure that the number of points has already been determined
		if(!determinedNumberOfPoints){
			System.out.println("Error: the number of points must be determined before reading the IFs as input.\nSystem.exit(0)");
			System.exit(0);
		}
		
		String line = null;
		int lineCounter = 0;
		int minIndex = 0;
		int maxIndex = 0;
		try{
			BufferedReader br = new BufferedReader(new FileReader(inputPath + filenameIFs));			
			while((line = br.readLine()) != null){
				String[] values = line.split("\t");
				if(values.length == 5){
					//  parse the values once that will be used to save time
					int firstIndex = Integer.parseInt(values[0]);
					int secondIndex = Integer.parseInt(values[1]);

					//  initialize the minimum index and the maximum index on the first round
					if(lineCounter == 0){
						minIndex = firstIndex;
						maxIndex = firstIndex;
					}

					//  check for a new minimum index
					if(firstIndex < minIndex){
						minIndex = firstIndex;
					}
					if(secondIndex < minIndex){
						minIndex = secondIndex;
					}

					//  check for a new maximum index
					if(firstIndex > maxIndex){
						maxIndex = firstIndex;
					}
					if(secondIndex > maxIndex){
						maxIndex = secondIndex;
					}
				}
				else{
					//  each row must have in order: index1    index2    IF    stdev    count
					System.out.println("Error: there are not five (5) entries in each row of the IFs file\n"
							+ "System.exit(0)");
					System.exit(0);
				}
				lineCounter++;
			}
			br.close();
		}catch (IOException e) {
			e.printStackTrace();
			System.out.println("Error: cannot find IFs file with filename \"" + filenameIFs + "\"\n" +
			"System.exit(0)");
			System.exit(0);
		}

		if(maxIndex != numberOfPoints){
			System.out.println("Error: the number of points is not the same in the IF values file (" + maxIndex 
					+ ") and the fragments file (" + numberOfPoints + ")\n"
					+ "System.exit(0)");
			System.exit(0);
		}


		IFvalues = new double[numberOfPoints][numberOfPoints];
		for(int i=0; i<numberOfPoints; i++){
			for(int j=0; j<numberOfPoints; j++){
				IFvalues[i][j] = -1;
			}
		}
		IFvariances = new double[numberOfPoints][numberOfPoints];
		for(int i=0; i<numberOfPoints; i++){
			for(int j=0; j<numberOfPoints; j++){
				IFvariances[i][j] = -1;
			}
		}

		line = null;
		try{
			BufferedReader br = new BufferedReader(new FileReader(inputPath + filenameIFs));			
			while((line = br.readLine()) != null){
				String values[] = line.split("\t");
				//  make sure that there are exactly five (5) entries in each row
				if(values.length == 5){
					int firstIndex = Integer.parseInt(values[0]) - minIndex;
					int secondIndex = Integer.parseInt(values[1]) - minIndex;
					
					//  exclude IF entries that do not have an associated standard deviation measure
					double stdev = Double.parseDouble(values[3]);
					if(stdev > 0){
						double variance = stdev * stdev;
						//  keep the IFvariances matrix symmetric
						IFvariances[firstIndex][secondIndex] = variance;
						IFvariances[secondIndex][firstIndex] = IFvariances[firstIndex][secondIndex];
						numberOfVarianceValues++;

						//  keep the IFvalues matrix symmetric
						IFvalues[firstIndex][secondIndex] = Double.parseDouble(values[2]);
						IFvalues[secondIndex][firstIndex] = IFvalues[firstIndex][secondIndex];
						averageIF += IFvalues[firstIndex][secondIndex];
						numberOfIFValues++;
					}									
				}
				else{
					System.out.println("Error: there are not five (5) entries in each row of the IFs file\n"
							+ "System.exit(0)");
					System.exit(0);
				}
			}
			br.close();
		}catch (IOException e) {
			e.printStackTrace();
			System.out.println("Error: cannot find IFs file with filename \"" + filenameIFs + "\"\n" +
			"System.exit(0)");
			System.exit(0);
		}
		System.out.println("Opened and read file: " + filenameIFs);

		
		//  calculate the average IF value
		if(numberOfIFValues == 0){
			System.out.println("Error: cannot divide by zero. numberOfIFValues is zero (0).\nSystem.exit(0)");
			System.exit(0);
		}
		else{
			averageIF = averageIF / numberOfIFValues;
		}
		if(numberOfIFValues != numberOfVarianceValues){
			System.out.println("Error: the number of IF entries (" + numberOfIFValues + ") does not equal the number of variance entries (" + numberOfVarianceValues + ")"
					+ "\nSystem.exit(0)");
			System.exit(0);
		}
		else{
//			System.out.println("the number of IF entries (" + numberOfIFValues + ") --- the number of variance entries (" + numberOfVarianceValues + ")");
		}
	}/*  end of getIFsFromFile method  */

	
	/*  start of initializeStructure method  */
	private void initializeStructure(BufferedWriter bwParameters) throws IOException{
		switch (initialStructureCase) {
		case 0:			//		INITIALIZE STRUCTURE FROM FILE
			structure = pdb.getStructure(inputPath + initialStructureFilename);
			bwParameters.write("Initializing the structure with a given structure from file \"" + initialStructureFilename + "\"");
			bwParameters.newLine();
			bwParameters.flush();
			break;
		case 1:			//		CUBE WITH SIDE LENGTH [0,0.01]
			structure = new double[numberOfPoints][3];

			Random random = new Random();
			
			//  initialize the structure in a cube with side length [0,0.01]
			for(int i=0; i<structure.length; i++){
				structure[i][0] = random.nextDouble() / 100.0;
				structure[i][1] = random.nextDouble() / 100.0;
				structure[i][2] = random.nextDouble() / 100.0;
			}
			bwParameters.write("Initializing the structure in a cube with side length [0,0.01]");
			bwParameters.newLine();
			bwParameters.flush();
			break;
		case 2:			//		CUBE WITH SIDE LENGTH = 10 * AVERAGE_DISTANCE
			structure = new double[numberOfPoints][3];

			random = new Random();
			
			double averageDistanceX10 = 10.0 / (Math.pow(averageIF, exponent));
			
			//  initialize the structure in a cube with side length = 10 * average distance of input dataset
			for(int i=0; i<structure.length; i++){
				structure[i][0] = random.nextDouble() * averageDistanceX10;
				structure[i][1] = random.nextDouble() * averageDistanceX10;
				structure[i][2] = random.nextDouble() * averageDistanceX10;
				
			}
			bwParameters.write("Initializing the structure in a cube with side length = " + averageDistanceX10);
			bwParameters.newLine();
			bwParameters.flush();
			break;
		case 3:			//		SPRIAL WITH SMALL RANDOM NOISE
			structure = new double[numberOfPoints][3];

			random = new Random();
			
			//  use counter1 for the x-coordinates 
			int counter1 = 1;
			//  use counte2 for the y-coordinates
			int counter2 = 2;
			//  use zStepSize for the z-coordinates
			double zStepSize = 5.0 / ((double)numberOfPoints);
			
			//  initialize the structure along a spiral with a bit of random noise
			for(int i=0; i<structure.length; i++){
				//  x-coordinate
				if(counter1 == 1 || counter1 == 2){
					structure[i][0] = 10 + random.nextDouble();
					counter1++;
				}
				else{
					structure[i][0] = 11 + random.nextDouble();
					counter1++;
				}
				//  reset the counter1
				if(counter1 == 5){
					counter1 = 1;
				}
				
				//  y-coordinate
				if(counter2 == 1 || counter2 == 2){
					structure[i][1] = 10 + random.nextDouble();
					counter2++;
				}
				else{
					structure[i][1] = 11 + random.nextDouble();
					counter2++;
				}
				//  reset the counter2
				if(counter2 == 5){
					counter2 = 1;
				}
				
				//  z-coordinate
				structure[i][2] = (i * zStepSize) + random.nextDouble();			
			}
			bwParameters.write("Initializing the structure in a spiral (with small random noise)");
			bwParameters.newLine();
			bwParameters.flush();
			break;
		default:		//  Need to specify an initial structure case
			System.out.println("Error: there was no initial structure case provided.\nSystem.exit(0)");
			System.exit(0);
			break;
		}
	}/*  end of initializeStructure method  */

	
	/*  start of outputGenomicPositionMap method  */
	private void outputGenomicPositionMap(){
		
		String outputFilename = outputPath + "genomicPositionMap_" + filenameIFs;

		try{						
			BufferedWriter bw = new BufferedWriter(new FileWriter(outputFilename));

			//  write the map
			for(int i=0; i<structure.length; i++){
				bw.write((i+1) + "\t" + pointGenomicLocations[i]);
				bw.newLine();
				bw.flush();				
			}
			bw.flush();
			bw.close();
		}catch (IOException e) {
			System.out.println("Error: could not output file: " + outputFilename);
			e.printStackTrace();
		}//end of try-catch block
		System.out.println("\tOutput file: " + outputFilename);
	}/*  end of outputGenomicPositionMap method  */

	
	/*  start of outputStructure method  */
	private void outputStructure(double structureIndexDouble, double probability, String name){		
		//  define the structure filename
		String outputFilename = new String();
		
		String structureIndexString = formatter.format(structureIndexDouble);
		
		if(name.compareTo("default") == 0){
			outputFilename += "structure_";
			
			//  determine the number of leading zeros
			if(structureIndexDouble < 10){
				outputFilename += "000";
			}
			else if(structureIndexDouble < 100){
				outputFilename += "00";
			}
			else if(structureIndexDouble < 1000){
				outputFilename += "0";
			}
			else{
				//  don't need to do anything
			}
			
			//  append the filetype extension
			outputFilename += structureIndexString + ".pdb";
		}
		else{
			//  use the name passed as a parameter to the method
			outputFilename += name;
		}
		
		
		//  create the header and write the structure to file
		String header = dateFormat.format(new Date()) + "\tlengthOfRun=\t" + formatter.format(runLength) + "\tstepSize=\t" + stepSize +
			"\tindex=\t" + structureIndexString + "\tln(probability)=\t" + probability;
		pdb.writeStructure(outputPath + outputFilename, structure, header);
	}/*  end of outputStructure method  */
	
	
	/*  start of runMCMC method  */
	private void runMCMC(BufferedWriter bwParameters) throws IOException{
		
		
		//  initialize the structure
		initializeStructure(bwParameters);
				

		//  output the genomic position map
		outputGenomicPositionMap();

		
		//  set the half step size
		double halfStepSize = stepSize / 2.0;

		
		//  create the random number generator
		Random random = new Random();

		
		//  do not impose any prior probabilities on the structures
		double probabilityOfStructure = 1.0;
		double probabilityOfIFData = 1.0;


		//  keep track of the best structure ever seen
		String bestStructureFilename = "best_structure_" + filenameIFs.substring(0, filenameIFs.lastIndexOf('.')) + ".pdb";
		bestStructure = new double[numberOfPoints][3];
		probabilityOfBestStructure = -1;
		indexOfBestStructure = -1;

		try{
			BufferedWriter bw = new BufferedWriter(new FileWriter(outputPath + "stats_" + filenameIFs));

			//  write the header line
			bw.write(outputFolder);
			bw.flush();
			bw.newLine();
			bw.write("command line args=\t");
			bw.flush();
			for(String arg: argsGlobal){
				bw.write(arg + "\t");
				bw.flush();
			}
			bw.newLine();

			
			//  calculate the probability of the initial structure
			//  ln(Pr[S|IF]) = ln(Pr[IF|S]) + ln(Pr[S]) - ln(Pr[IF])
			//	ln(Pr[IF|S]) = SUM(i,j)[- ((IF(i,j) - (1/d(i,j)))^2 / variance)]
			
			double probabilityCurrent = 0;
			
			for(int i=0; i<numberOfPoints-1; i++){
				for(int j=i+1; j<numberOfPoints; j++){
					if(IFvalues[i][j] != -1.0){

						//  calculate the distance in the structure between point i and point j
						double distanceTerm1 = structure[j][0] - structure[i][0];
						distanceTerm1 *= distanceTerm1;
						double distanceTerm2 = structure[j][1] - structure[i][1];
						distanceTerm2 *= distanceTerm2;
						double distanceTerm3 = structure[j][2] - structure[i][2];
						distanceTerm3 *= distanceTerm3;
						
						double distance = distanceTerm1 + distanceTerm2 + distanceTerm3;
						distance = Math.sqrt(distance);
						
						
						
						
						//  before switch (2011-March-18):
//						double IFmodel = 1.0 / distance;
						//  after switch (2011-March-18):
						double IFmodel = 1.0 / Math.pow(distance, exponent);

						
						
						
						//  calculate the two terms and add the product to the sum
						double term = IFvalues[i][j] - IFmodel;	//  the numerator
						term *= term;
						term /= IFvariances[i][j];		//  divide the numerator by the denominator
						term += Math.log(2.0 * Math.PI * IFvariances[i][j]);

						probabilityCurrent += term;
					}
				}
			}
			probabilityCurrent *= -0.5;
			
			//  output the initial structure
			outputStructure(0, probabilityCurrent, "structure_0000.pdb");
			bw.write(0 + "\t" + probabilityCurrent);
			bw.newLine();
			bw.flush();

			
			//  set the initial structure as the best structure seen
			
			probabilityOfBestStructure = probabilityCurrent;
			indexOfBestStructure = 0;
			for(int i=0; i<structure.length; i++){
				bestStructure[i][0] = structure[i][0];
				bestStructure[i][1] = structure[i][1];
				bestStructure[i][2] = structure[i][2];
			}
			

			//  use a separate counter so the structures that are output have consecutive numbers in their filenames
			int structureFilenameCounter = 1;
			
			
			//  run the Markov chain Monte Carlo simulation using the Metropolis-Hastings algorithm
			for(long counter=1; counter<=runLength; counter++){

				//  choose a point to move at random
				int pointToMove = random.nextInt(numberOfPoints);
				double currentX = structure[pointToMove][0];
				double currentY = structure[pointToMove][1];
				double currentZ = structure[pointToMove][2];

				
				//  move the point a random distance
				double newX = currentX + (random.nextDouble() * stepSize) - halfStepSize;
				double newY = currentY + (random.nextDouble() * stepSize) - halfStepSize;
				double newZ = currentZ + (random.nextDouble() * stepSize) - halfStepSize;
				

				//  calculate the probability of the new sampled structure
				//  ln(Pr[S|IF]) = ln(Pr[IF|S]) + ln(Pr[S]) - ln(Pr[IF])
				//	ln(Pr[IF|S]) = SUM(i,j)[- ((IF(i,j) - (1/d(i,j)))^2 / variance)]
				double probabilityNew = probabilityCurrent;
								
				for(int j=0; j<numberOfPoints; j++){
					if(j != pointToMove){
						if(IFvalues[pointToMove][j] != -1.0){
							double xj = structure[j][0];
							double yj = structure[j][1];
							double zj = structure[j][2];
							
							//  calculate the distance in the current structure between the point to move (current location) and point j
							double distanceCurrentTerm1 = xj - currentX;
							distanceCurrentTerm1 *= distanceCurrentTerm1;
							double distanceCurrentTerm2 = yj - currentY;
							distanceCurrentTerm2 *= distanceCurrentTerm2;
							double distanceCurrentTerm3 = zj - currentZ;
							distanceCurrentTerm3 *= distanceCurrentTerm3;

							double distanceCurrent = distanceCurrentTerm1 + distanceCurrentTerm2 + distanceCurrentTerm3;
							distanceCurrent = Math.sqrt(distanceCurrent);
							
							//  make sure not to divide by zero
							if(distanceCurrent <= 0){
								System.out.println("Error: the distance between the points was not properly determined\n"
										+ "System.exit(0)");
								System.exit(0);
							}
							
							
							
							
							
							
							//  before switch (2011-March-18):
//							double IFmodelCurrent = 1.0 / distanceCurrent;
							//  after switch (2011-March-18):
							double IFmodelCurrent = 1.0 / Math.pow(distanceCurrent, exponent);

							
							
							
							
							
							//  calculate the distance in the new structure between the point to move (new location) and point j
							double distanceNewTerm1 = xj - newX;
							distanceNewTerm1 *= distanceNewTerm1;
							double distanceNewTerm2 = yj - newY;
							distanceNewTerm2 *= distanceNewTerm2;
							double distanceNewTerm3 = zj - newZ;
							distanceNewTerm3 *= distanceNewTerm3;

							double distanceNew = distanceNewTerm1 + distanceNewTerm2 + distanceNewTerm3;
							distanceNew = Math.sqrt(distanceNew);
							
							//  make sure not to divide by zero
							if(distanceNew <= 0){
								System.out.println("Error: the distance between the points was not properly determined\n"
										+ "System.exit(0)");
								System.exit(0);
							}
							
							
							
							
							
							
							//  before switch (2011-March-18):
//							double IFmodelNew = 1.0 / distanceNew;
							//  after switch (2011-March-18):
							double IFmodelNew = 1.0 / Math.pow(distanceNew, exponent);
							
							
							
							
							
							
							//  optimized fast method (used up until 18-March-2011
							//  update the probability of the new structure (only change the terms that contain the point to move)
//							double factor1 = 1.0 / (2 * IFvariances[pointToMove][j]);
//							double factor2 = IFmodelCurrent - IFmodelNew;
//							double factor3 = IFmodelNew + IFmodelCurrent - (2 * IFvalues[pointToMove][j]);
//							double factorUpdate = factor1 * factor2 * factor3;
//							
//							probabilityNew += factorUpdate;
							
							
							//  slower method
							//  update the probability of the new structure (only change the terms that contain the point to move)							
							//  first subtract the value of the term corresponding to the current location of the point
//							probabilityNew -= (-1.0 / (2.0 * IFvariances[pointToMove][j])) * Math.pow(IFmodelCurrent - IFvalues[pointToMove][j], 2);
							//  second add the value of the term corresponding to the new location of the point
//							probabilityNew += (-1.0 / (2.0 * IFvariances[pointToMove][j])) * Math.pow(IFmodelNew - IFvalues[pointToMove][j], 2);
							
							
							//  non-optimized method implemented to allow for different exponents in the conversion of distance to IF
							//  used starting 18-March-2011
							//  note: this is the "slower method" described above
							double varianceTerm = -1.0 / (2.0 * IFvariances[pointToMove][j]);
							double IFvalueExperimental = IFvalues[pointToMove][j];
							
							//  first calculate the term corresponding to the old model
							double IFdifferenceOldSquared = IFmodelCurrent - IFvalueExperimental;
							IFdifferenceOldSquared *= IFdifferenceOldSquared;

							//  second calculate the term corresponding to the new model
							double IFdifferenceNewSquared = IFmodelNew - IFvalueExperimental;
							IFdifferenceNewSquared *= IFdifferenceNewSquared;
							
							//  third calculate the difference between the old model and the new model
							double IFnetDifference = IFdifferenceNewSquared - IFdifferenceOldSquared;
							
							//  calculate the probability of the new structure
							probabilityNew += varianceTerm * IFnetDifference;
						}
					}
				}


				//  accept or reject the new structure using the acceptance function
				double lnProbabilityRatio = probabilityNew - probabilityCurrent;
				double probabilityRatio = Math.pow(Math.E, lnProbabilityRatio);
				double alpha = 1.0;
				if(probabilityRatio < alpha){		//  if the ratio is less than one (1), then the new structure has a lower probability than the current structure
					alpha = probabilityRatio;
				}

				if(alpha >= random.nextDouble()){				//  accept the new structure if alpha is bigger than a uniformly chosen random variable [0,1[
					//  update the probability of the structure					
					probabilityCurrent = probabilityNew;
					
					
					//  keep track of the number of changes accepted during the MCMC run
					numberOfChangesAccepted++;


					//  accept the new structure by updating the coordinates of the point to move
					structure[pointToMove][0] = newX;
					structure[pointToMove][1] = newY;
					structure[pointToMove][2] = newZ;

					
					//  check if the new structure is the best structure ever seen
					if(probabilityCurrent > probabilityOfBestStructure){
						probabilityOfBestStructure = probabilityCurrent;
						indexOfBestStructure = counter;

						//  update the coordinates of the best structure
						for(int i=0; i<structure.length; i++){
							bestStructure[i][0] = structure[i][0];
							bestStructure[i][1] = structure[i][1];
							bestStructure[i][2] = structure[i][2];
						}
					}
				}

				//  output the structure and statistics to file
				if((counter % (runLength / numberOfStructuresToOutput)) == 0){				
					outputStructure(structureFilenameCounter, probabilityCurrent, "default");
					bw.write(structureFilenameCounter + "\t" + probabilityCurrent);
					bw.flush();
					bw.newLine();
					structureFilenameCounter++;
				
					//  output and overwrite the best structure seen so far
					double[][] tempStructure = structure;
					structure = bestStructure;				//  note: the outputStructure method will output the structure pointed to by the structure global variable
					outputStructure(indexOfBestStructure, probabilityOfBestStructure, bestStructureFilename);
					structure = tempStructure;					
				}
			}//end of Markov chain Monte Carlo loop
			bw.close();
		}catch (IOException e) {
			System.out.println("Error: could not output file: " + outputPath + "stats_" + filenameIFs);
			e.printStackTrace();
		}//end of try-catch block
		System.out.println("\tOutput file: " + outputPath + "stats_" + filenameIFs);
		
		
		//  output the best structure seen during the run
		structure = bestStructure;				//  note: the outputStructure method will output the structure pointed to by the structure global variable
		outputStructure(indexOfBestStructure, probabilityOfBestStructure, bestStructureFilename);
		bwParameters.write("bestStructureProbability=\t" + probabilityOfBestStructure);
		bwParameters.newLine();
		bwParameters.flush();
	}/*  end of runMCMC method  */


	/*  start of run method  */
	public void run(String[] args) {
		//  make the command line arguments accessible by all methods as a global variable
		argsGlobal = args;
		
		switch (args.length) {
		case 2:
			filenameIFs = args[0];
			filenameFragments = args[1];
			break;
		case 3:
			filenameIFs = args[0];
			filenameFragments = args[1];
			runLength = Long.parseLong(args[2]);
			break;
		case 4:
			filenameIFs = args[0];
			filenameFragments = args[1];
			runLength = Long.parseLong(args[2]);
			numberOfStructuresToOutput = Double.parseDouble(args[3]);
			break;
		case 5:
			filenameIFs = args[0];
			filenameFragments = args[1];
			runLength = Long.parseLong(args[2]);
			numberOfStructuresToOutput = Double.parseDouble(args[3]);
			stepSize = Double.parseDouble(args[4]);
			break;
		case 6:
			filenameIFs = args[0];
			filenameFragments = args[1];
			runLength = Long.parseLong(args[2]);
			numberOfStructuresToOutput = Double.parseDouble(args[3]);
			stepSize = Double.parseDouble(args[4]);
			exponent = Double.parseDouble(args[5]);
			break;
		case 7:
			filenameIFs = args[0];
			filenameFragments = args[1];
			runLength = Long.parseLong(args[2]);
			numberOfStructuresToOutput = Double.parseDouble(args[3]);
			stepSize = Double.parseDouble(args[4]);
			exponent = Double.parseDouble(args[5]);
			initialStructureCase = Integer.parseInt(args[6]);
			break;
		case 8:
			filenameIFs = args[0];
			filenameFragments = args[1];
			runLength = Long.parseLong(args[2]);
			numberOfStructuresToOutput = Double.parseDouble(args[3]);
			stepSize = Double.parseDouble(args[4]);
			exponent = Double.parseDouble(args[5]);
			initialStructureCase = Integer.parseInt(args[6]);
			initialStructureFilename = args[7];
			break;
		default:
			System.out.println("Error:\tplease use one of the following patterns (without the quotes)\n" +
					"\t\"java MCMC5CMain IFsFilename.txt fragmentsFilename.txt\"\n" +
					"\t\"java MCMC5CMain IFsFilename.txt fragmentsFilename.txt runLength\"\n" +
					"\t\"java MCMC5CMain IFsFilename.txt fragmentsFilename.txt runLength numberOfStructuresToOutput\"\n" +
					"\t\"java MCMC5CMain IFsFilename.txt fragmentsFilename.txt runLength numberOfStructuresToOutput stepSize\"\n" +
					"\t\"java MCMC5CMain IFsFilename.txt fragmentsFilename.txt runLength numberOfStructuresToOutput stepSize exponent\"\n" +
					"\t\"java MCMC5CMain IFsFilename.txt fragmentsFilename.txt runLength numberOfStructuresToOutput stepSize exponent initialStructureCase\"\n" +
					"\t\"java MCMC5CMain IFsFilename.txt fragmentsFilename.txt runLength numberOfStructuresToOutput stepSize exponent initialStructureCase initialStructureFilename\"\n" +
				"Notes:\tthe default run length is " + formatter.format(runLength) + "\n" +
					"\tthe default number of structures to output is " + (int)numberOfStructuresToOutput + "\n" +
					"\tthe default step size is " + stepSize + "\n" +
					"\tthe default exponent is " + exponent + "\n" +
					"\tthe default initial structure case is " + initialStructureCase + "\n" +
					"\tthe default initial structure filename is \"" + initialStructureFilename + "\"\n" +
				"System.exit(0)");
			System.exit(0);
			break;
		}

		
		//  create the output folder
		new File(outputFolder).mkdir();
		
		
		//  get the number of points (can pass the fragments file as parameter to count the number of points from this file)
		//  the number of points MUST be determined before the fragments or the IFs are read as input
		getNumberOfPoints(filenameFragments);

		
		//  get the fragment genomic locations
		getFragmentsFromFile();


		//  get the IF values (must be done after getting the fragment genomic locations)
		getIFsFromFile();


		//  calculate how long running the MCMC takes (time point before run)
		long timeBefore = System.currentTimeMillis();


		//  output the parameters to file
		try{
			String parameterFilename = "parameters_" + filenameIFs;
			BufferedWriter bwParameters = new BufferedWriter(new FileWriter(outputPath + parameterFilename));
			
			bwParameters.write(outputPath);
			bwParameters.newLine();
			bwParameters.flush();
			bwParameters.write("filenameIFs=\t" + filenameIFs);
			bwParameters.newLine();
			bwParameters.flush();
			bwParameters.write("filenameFragments=\t" + filenameFragments);
			bwParameters.newLine();
			bwParameters.flush();
			bwParameters.write("runLength=\t" + runLength);
			bwParameters.newLine();
			bwParameters.flush();
			bwParameters.write("numberOfStructuresToOutput=\t" + (int)numberOfStructuresToOutput);
			bwParameters.newLine();
			bwParameters.flush();
			bwParameters.write("stepSize=\t" + stepSize);
			bwParameters.newLine();
			bwParameters.flush();
			bwParameters.write("exponent=\t" + exponent);
			bwParameters.newLine();
			bwParameters.flush();
			bwParameters.write("command line parameters:");
			if(args.length == 0){
				bwParameters.write("\tnone");
			}
			else{
				for(String arg: args){
					bwParameters.write("\t" + arg);
				}
			}
			bwParameters.newLine();
			bwParameters.flush();
			
			
			
			
			//  run the MCMC simulation
			runMCMC(bwParameters);
			
			
			
			
			//  output the time taken by the MCMC run
			long timeTaken = System.currentTimeMillis() - timeBefore;
			timeTaken = timeTaken / (long)(1000.0);			//	convert from milliseconds to seconds		
			if(timeTaken > 3600){		//	(1hr = 3600 sec)
				int numberOfHours = (int)timeTaken / (int)3600;
				long leftOver = timeTaken % 3600;
				int numberOfMinutes = (int)leftOver / (int)60;
				leftOver = leftOver % 60;
				int numberOfSeconds = (int)leftOver;
				bwParameters.write("Amount of time taken = " + numberOfHours + "hr" + numberOfMinutes + "min" + numberOfSeconds + "sec");
				bwParameters.newLine();
				bwParameters.flush();
				//				System.out.println("Amount of time taken = " + numberOfHours + "hr" + numberOfMinutes + "min" + numberOfSeconds + "sec");				
			}
			else if(timeTaken > 60){		//	(1min = 60 sec)
				int numberOfMinutes = (int)timeTaken / (int)60;
				long leftOver = timeTaken % 60;
				int numberOfSeconds = (int)leftOver;
				bwParameters.write("Amount of time taken = " + numberOfMinutes + "min" + numberOfSeconds + "sec");
				bwParameters.newLine();
				bwParameters.flush();
//				System.out.println("Amount of time taken = " + numberOfMinutes + "min" + numberOfSeconds + "sec");
			}
			else{
				int numberOfSeconds = (int)timeTaken;
				bwParameters.write("Amount of time taken = " + numberOfSeconds + "sec");
				bwParameters.newLine();
				bwParameters.flush();
//				System.out.println("Amount of time taken = " + numberOfSeconds + "sec");
			}
			
			
			//  output the number of proposed changes accepted during the run
			bwParameters.write("Number of proposed changes accepted: " + numberOfChangesAccepted +
					" [" + (int)(numberOfChangesAccepted * 100 / runLength) + "%]");
			bwParameters.newLine();
			bwParameters.flush();
//			System.out.println("Number of proposed changes accepted: " + numberOfChangesAccepted +
//					" [" + (int)(numberOfChangesAccepted * 100 / runLength) + "%]");
			
			
			bwParameters.close();
			
			System.out.println("\tOutput file: " + outputPath + parameterFilename);
		}catch (IOException e) {
			e.printStackTrace();
		}//end of try-catch block
		
		System.out.println("Done!");
	}/*  end of run method  */
}/*  end of MCMC5C class  */